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In this paper the internal structure of a neutron star is shown to be inferrable from its 
gravitational-wave spectrum. Iteratively applying the inverse scheme of the scaled coordinate log- 
arithmic perturbation method for neutron stars proposed by Tsui and Leung [Astrophys. J. 631, 
495 (2005)], we are able to determine the mass, the radius and the mass distribution of a star from 
its quasi-normal mode frequencies of stellar pulsation. In addition, accurate equation of state of 
nuclear matter can be obtained from such inversion scheme. Explicit formulas for the case of axial 
ui-mode oscillation are derived here and numerical results for neutron stars characterized by different 
\Q • equations of state are shown. 
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O . I. INTRODUCTION 

O; 

As remnants of supernova explosions, neutron stars comprised of matters with subnuclear and supranuclear densities 
£SJ ' have long been a major subject of intense interest for researchers in astrophysics, nuclear physics and particle physics 
since the pioneering studies of Oppenheimer and Volkoff [l| and Tolman |2(. Their internal structure is likely to 
provide a direct test for theories of nuclear matter, quark matter and high energy physics (see e.g. [J and 
references therein). For example, through measurement of the Eddington flux and the redshift in the X-ray burst, the 
mass and the radius of the neutron star EXO 0748-676 were simultaneous determined, and relevant data were used to 
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rule out soft equations of state (EOS) for nuclear matter [EEl- Likewise, a few astronomical observations revealing 
evidences for the existence of the strange quark star (SOS), a variant of the neutron star (see e.g. [EIEE3)j to some 
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extent resuscitate the theory of quark matter [Til [ 

On the other hand, the huge gravitational field accompanying a dense neutron star also attracts the attention 
of relativists. In particular, undulating neutron stars are expected to be promising sources of gravitational waves. 
For instance, when gravitational-wave detectors of various designs are in full swing in one or two decades (see e.g. 
[TH Il6| and references therein), the frequency of detection of g ravitational waves emitted in the mergers of binary 
5^ \ neutron stars could be as high as several hundreds per year |!5lll7l ll8 | . Besides, gravitational waves emitted during 
&)[)■ asymmetric stellar core collapses leading to the formation of neutron stars could also be detectable 0, |2(j ■ It is then 
likely that neutron stars could be surveyed in the gravitational wave channel. More importantly, gravitational waves 
can easily penetrate the core of a neutron star without being absorbed, and might carry useful information about the 
deep interior of the star. As the behavior of nuclear matter at high densities is not yet clearly known, the internal 
structure of neutron stars is masked and often gives rise to many debatable issues. Hence, a lot of in-depth studies have 
been sparked off by the possibility of unveiling the internal structure of neutron stars from their gravitational-wave 
signals. Detailed examination of such signals, which has been coined "gravitational- wave asteroseismology" , were 
carried out by various groups of researchers (see e.g. 0, HE HE 0> HE HE HE HE H3) to study effects of different 
physical processes, including superfluidity, occurrence of quark matter, and phase transitions. 

Gravitational waves emitted from a neutron star are commonly analyzed in terms of quasi-normal modes (QNMs) 
characterized by complex eigenfrequencies w = cu r — iu>i, with u> r and u>i measuring respectively the rates of pulsation 
and damping |30l l31| . To infer the physical parameters of a neutron star from its QNM frequencies, Andersson and 
Kokkotas made use of the universal behavior in the QNMs of the fundamental fluid /-mode and the first polar to-mode 
to estimate the radius and the mass of the star Besides, the most probable EOS among all the known models of 
stellar matters could also be identified from the frequency of the leading p-mode (or axial w-mode) oscillation p^.l23|. 

In a series of recent papers [IE Ull, we have studied the physical or igin of the universality in the pulsation frequencies 
of the /-mode and if- mode stellar pulsations discovered in Refs. ^2 an d HE- an d developed a scaled-coordinate 
logarithmic perturbative theory (SCLPT) for neutron stars to locate axial w-mode oscillations. The main objective 
of the current paper is then to work out a feasible scheme to infer the internal structure of a neutron star from its 
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gravitational wave spectrum. Assuming that the frequencies and the damping rates of a few w-mode oscillations 
of a neutron star can be identified from gravitational-wave observation, we iteratively invert the SCLPT scheme to 
determine the mass M, the radius R, the mass density p(r) at radius r, and the mass distribution function m(r) (i.e. 
the total mass inside radius r) of the star. In turn, accurate equation of state of nuclear matter can also be obtained 
numerically from the mass distribution function. While the idea underlying the inversion scheme reported here is 
completely generic, we will consider and use the spectrum of axial w-mode oscillations to illustrate our method. As a 
remark, we note here that Lindblom has proposed a scheme to determine, up to certain accuracy, EOS of the nuclear 
matter of a neutron star from its mass and radius [34j . However, the EOS considered in Lindblom's scheme is limited 
to barotropic type. In contrast, our method directly exploit the frequencies of several QNMs to determine the radius, 
the mass and other useful parameters of the star. Moreover, the EOS obtained from our inversion scheme is free from 
constraints. We will show in the following discussion that even the exotic EOS of quark matter can be inferred from 
the gravitational wave signals of a SQS. 

The organization of our paper is as follows. In Sect. II, we review the universality existing in the QNM frequencies 
(including the real and the imaginary parts) of neutron stars described by different EOS and discuss how such behavior 
can be explained by the Tolman VII model (TVIIM) of stars Q. In Sect. Ill we approximate the mass distribution 
of a realistic neutron star by that of TVIIM plus some small correction terms to be determined from the inversion 
scheme. Section IV outlines the spirit of the inversion scheme while explicit formulas for the scheme are derived in 
Sect. V. In Sect. VI we apply our method to invert the EOSs of several realistic stars and numerical results are shown 
there. Section VII studies the feasibility of our method and then we conclude our paper in Sect. VIII. The present 
paper is the expanded version of a recent Letter [3f|. Unless otherwise stated, geometric units are used in this paper. 



II. UNIVERSALITY AND TVIIM 



First of all, we review the universality in QNMs of gravitational waves and state its physical origin, thus introducing 
the starting point of our inversion scheme — the TVIIM. It has recently been observed that the QNM frequencies of 
the leading u>-mode and /-mode oscillations of non-rotating neutron stars approximately satisfy the following scaling 
law [2II23I: 



Mlo 



(2.1) 



Here a, 6, and c are complex constants determined from curve fitting. In Fig. ^ we plot Mu> r and Mu>i against the 
compactness C = Mj R for the leading axial ui-mode of non-rotating neutron stars characterized by different EOSs, 
including models A and C proposed by Pandharipande |3(| |13 , models AU and UT proposed by Wiringa, Fiks and 
Fabrocini (3^, models APR1 and APR2 proposed by Akmal, Pandharipande and Ravenhall One can easily see 
that both Mu> r and MtOi are, in agreement with i|2.1J) . well approximated by quadratic function of C. 

The universality observed in the spectra of realistic neutron stars comes as a surprise because, as shown in Fig. 
the pressure P at a given density p varies greatly for the EOSs considered there. To understand the physical origin 
underlying the universality, we showed in Ref. J33 that the scaled complex eigenfrequencies Mlo of axial w-mode, 
polar w-mode and /-mode oscillations, to a good approximation, depend only on the compactness M/R. In addition, 
we also discovered that the TVIIM, which was first proposed by Tolman in 1939 |2| and has a mass distribution 
function given by: 



m c (r) = M 
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can mimic the mass distribution of most neutron stars. As clearly shown in Fig. EI where the normalized mass function 
m(r)/M is plotted against the scaled radius r/R for stars with a compactness M/R = 0.28 and described by different 
EOS as mentioned above, the mass distribution function m c (r) is able to approximate the mass distribution inside 
neutron stars with different EOSs. Moreover, the scaled frequency Mlo of QNMs of TVIIM is close to those of realistic 
neutron stars (see Fig.^l and display similar scaling behavior mentioned above. It is obvious that TVIIM does provide 
a good approximation to stars with varying EOSs, and the universality shown in Fig. ^ can be captured by the best 
quadratic fit to the QNMs of TVIIM, with a = -4.4 - 6.3i, 6 = 3.1 + 1.9i, and c = -0.072 + 0.098L A distinct feature 
of TVIIM is that the scaled complex eigenfrequencies Mlo of its QNMs depend only on the compactness of the star 
[3^|. Hence, the observed universality in the QNMs (w-mode and /-mode) is expected to be closely related to TVIIM 
and was indeed verified using SCLPT [33| . 
The density distribution of TVIIM, given by 



P(r) = Pa 
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with the central density po — 15M/(8ttR 3 ), is so simple that the metric coefficients e^ r ' and e xtyr \ defined by the line 
element ds: 



can be obtained analytically: 



Here £ = r/R, and 



ds 2 = - e "(r) dt 2 + e A(r) dr 2 + r 2 (d(? 2 + 

e~ x = 1 -C£ 2 (5-3£ 2 ), 
e v = (1 - 5C/3)cos 2 0. 

4> = (w-l- w)/2 + 0i, 



sin 



W = l0E 



3C 



0i = 0(£ 2 = 1) = arctan y ^ _ , 
Wl = w(£ 2 = 1) . 

Besides, the pressure can also be found exactly from the Tolman-Oppenheimer-Volkoff (TOV) equations 0,0 : 
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V3Ce~ A tan0- -(5- 3£ 2 ) 



By virtue of (|2.3fl . we have: 

e" A = 1 - C ( 1 
The EOS of TVIIM is consequently given by: 
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It is obvious that both m c (r) and the EOS for TVIIM depend on M and R (or any two independent parameters 
derivable from them). Thus, once these two parameters are known, the mass distribution inside the star and the EOS 
of the stellar matter are determined automatically. In the following discussion, we will make use of this salient feature 
of TVIIM as the foundation of the inversion scheme. 



III. APPROXIMATE MASS DISTRIBUTION 

Based on the the universality and the TVIIM discussed above, a feasible scheme is proposed here to infer the 
internal structure of a neutron star from its gravitational wave spectra. As the first step of the scheme, we consider 



the frequency of the leading (i.e. the least-damped) axial u>-mode for TVIIM, uj^ — uj y {"J — iu))^' , and plot the ratio 
Q = uj^ /u>$ against with the compactness in Fig. 0J We see that lu[ c J /w[f is in fact a monotonically function of 
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M/R. Hence, once cj^ is known, M/R can be obtained from Fig. 0] and in turn Midf 1 can be read from Fig. ^ It 

(c) (c) 

is then straightforward to find M and R. Likewise, M and R can be inferred from ui 2 and oj^ , the frequencies of 
the second and the third leading modes, as an d w 3rV w 3? are a l so monotonic functions of the compactness 

C (see Fig.0J). However, the frequency u^' is preferred in the subsequent discussion because it has the most sensitive 
dependence on C. 

As TVIIM provides a benchmark for other realistic stars, we expect that the eigenfrequency of the leading axial 
w-mode emitted from a realistic neuron star, u}\ = u>i r — ioJn, is close to that of a TVIIM star with the same mass M 
and the same radius R. Therefore, we could go through the procedure outlined in the last paragraph with u>\° replaced 
by ll>i to obtain an estimate of the mass and the radius of the star in consideration. This first step in our scheme is 
analogous to the method proposed by Andersson and Kokkotas 22]. However, we will consider the frequencies of the 
non-leading modes and show that they can lead to much improved estimates of the mass, the radius and the EOS as 
well. 

The mass distribution of a realistic neutron star is of course different from that of a TVIIM star. However, the 
difference is small. Considering the asymptotic behavior of the density around r = (see, e.g. |39j), we expand m(r) 
as: 

p-2 

m(r, { Pk }, M, R) « m (r, M, R) + J2 W m j( r > M > R )> 

where 

mi(r ,M, J R) = ^(^) 2j (r 2 - J R 2 ) 2 , (3.2) 

p > 2, and iij (j = 1, 2, . . . ,p — 2) are p — 2 adjustable parameters to provide the best fit to m(r, {/J-k}, M, R). As the 
mass M and the radius R of the neutron star are not yet exactly known, they are also considered as free parameters. 
There are consequently p free parameters to be determined from the QNM frequencies. Besides, as the zeroth order 
approximation to m(r), mo is well approximated by m c (r), which will be chosen to be the initial guess for mo. 

Generally speaking, we can approximate the configuration of a neutron star by choosing a set of fij, M and R 
which "best fit" the exact one. Before working out the optimization process that determines these p parameters, we 
briefly mention the physical significance of m,j(r, M, R). It is clear that rrij{r — R,M,R) = for all j. Hence, rrij 
does not affect the total mass of the neutron star while modifying the internal mass distribution. Furthermore, it 
can be shown that each rrij has only one maximum value and the position of such maximum gets closer to the stellar 
surface as j increases. This implies that rrij with small j has significant contribution to the core of the neutron star, 
while rrij with large j mainly affects the the mass distribution near the stellar surface. In addition, the central density 
of the neutron star depends only on mo and mi. As a consequence, as far as the high density regime of the EOS is 
concerned, rrij with smaller j is likely to be more important. Lastly, the density distribution can be easily derived 
from the mass distribution, 

p-2 

P (r, n k ,M, R) = p {r, M, R) + ^ fi jPj (r, M, R), (3.3) 

j=i 

where 

15M 



po{r,M, R) 



Pj (r, M,R) = —=(-) (r 2 - R 2 ) [(5 + 2j)r 2 (1 + 2j)R 2 } . (3.5) 
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IV. THE INVERSION SCHEME 
A. Optimal model 



Suppose that the mass and the radius of a star are known with certain accuracy from the the frequency of its 
leading mode, which are equal to Mg and Rq, respectively. The exact values of these quantities, M = Mq + M\ and 
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R = Ro + R\ , and the mass distribution inside the star can be determined from the frequencies of QNMs of higher 
orders. Suppose n (n > 1) QNM frequencies of a neutron star are known and we are going to search for the optimal 
stellar model that yields the minimal deviation from the measured QNM frequencies. As discussed above, there are 
p independent small parameters in the model, namely fij (j — 1, 2, ... ,p — 2), /i p -i = M/R — Mq/Rq = , and 

fip = Mi . If all of these free parameters equal zero, the frequencies are denoted by u> q ^ . 

Based on the first order result of SCLPT [H], the corresponding first order change in the g-th (q = 1,2, ... ,n) 
QNM frequencies io q is 



^=I>^=I>^, (4-1) 

3=1 ^ 3 = 1 



which is close to the exact change: 



A q =u; q -J°K (4.2) 



Explicit expressions for c q j = dLu q /dfj,j will be derived in the following section of this paper using SCLPT. By 
minimizing the square-deviation 



x 2 = Ei*,-^ 5 

9=1 

= A i^i ~ A iY1 frCtf ~ A qYl MiC«i + ^43 1 » ( 4 - 3 ) 



9=1 

p P p 



q—1 \ j — 1 i—1 i—l j — 1 



where the overbar indicates the complex conjugate of a variable, one can obtain a set of /ij and hence an optimal 
approximate model for the observed neutron star. The necessary condition for \ 2 to achieve a minimum is obviously: 

= 0, (4.4) 

for all k, leading to n linear equations: 

n p n 

/I {AqCqk + \c qk ) = 22fJ,j 2J (CqkCqj + C qk C q0 ) . (4.5) 
9=1 3 = 1 9=1 

This system of equations can be rewritten in the following matrix representation: 

b = A(2, (4.6) 
where A is a p x p square matrix, fl is the column vector of {fJ,j}, and 

n 

b k = ^ (A q c qk + A q c qk ) , (4.7) 

9=1 

n 

A k3 = CqkCqj + C qk C qj ) , (4.8) 

9=1 

for k, j — 1, 2, 3, . . . ,p. Consequently, we can compute the best fij by solving the above matrix equation. 

B. Iterative scheme 

With the p parameters {^} obtained through the scheme sketched above, we can construct a new model star 
by incorporating these mass corrections into mo(r) and compute a new set of QNM frequencies for such a mass 
distribution. Again we can compare these QNM frequencies with the observed QNM frequencies uj q , and minimize 
the error % 2 by solving l|4.6|l to yield a new set of {/ij}. Through this iterative scheme, one can generate convergent 
sequences for the mass, the radius and the mass distribution. 
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C. Inversion of EOS 



The equilibrium configuration of a neutron star is given by the TOV equations P, : 

^ = 4nr 2 p, (4.9) 
dr 

dv 2m + 8Trr 3 P 

-r = , „ s . (4-10) 
dr r(r — 2m) 

Here p and P are evaluated at a point (r,0, Once m(r) is known, the density p(r) can be found from Q4.9JI . 
Combining (|4.10() and l|4.11|l to yield a first-order differential equation for P(r): 

d-P _^ p^m + Anr 3 P (4 12) 

dr r(r — 2m) ' 

one can obtain P(r) by solving (|4.12|l with the boundary condition P(r = R) = 0. Hence, the EOS can be straight- 
forwardly inferred by comparing the two functions p(r) and P(r). 

V. SCLPT 
A. Scaled Regge- Wheeler equation 

The starting point of the present inversion scheme is (|4.1() , which expresses the first-order shift in QNM frequencies 
as a linear sum over various perturbations {//j}. To evaluate the coefficients {c q j} in H4.ll) . we have to invoke the 
theory of first-order SCLPT to find how various perturbations {p,j} affect the QNM frequency [8^ |. 

First of all, we introduce the following scaled variables: 

(5.1) 
(5.2) 

(5.3) 

(5.4) 
(5.5) 
(5.6) 
(5.7) 
(5.8) 





r 


r 


M ' 




r* 


f* 


M ' 


m(f) 


m(r) 


M ' 


P(f) 


= M 2 P(r), 


p(f) 


= M 2 p(r), 


u(f) 


= v{r), 


Q 


= Mco, 


V{h) 


= M 2 V(r*) 



Here the tortoise coordinate r* is related to the circumferential radius r by 

-" +A)/2 dr, (5.9) 
and outside the star, the tortoise radial coordinate reduces to 



o 



r„=r + 2Mln(^-l)+C, (5.10) 

where C 1 is a constant that can be obtained by matching (|5.9() with (|5.10() at r = R. Besides, we also let R = R/M 
and P* = R*/M, where P* = r*(r = R) is the tortoise radius of the star. In terms of these scaled variables, the 
neutron star Regge- Wheeler (NSRW) equation j3^ governing the evolution of axial w-mode pulsation can be rewritten 
as [33|: 

d 2 1 ~ 

— +Q 2 -v(n) V(n) = o, (5.ii) 

art J 
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where for r < R the scaled Regge- Wheeler potential takes the form: 



= + l)f + 47rf J ( / 5 - P) - 6m(r)], 



(5.12) 



and for r > R, 



1(1 + 1) 6 



(5.13) 



The scaled NSRW equation l|5.11|l is obviously more amenable to perturbative expansion since outside the star V(f) 
depends only on f. Hence, changes in the stellar mass do not directly affect the potential outside the star. Therefore, 
perturbation arising from the change in the stellar mass, which also leads to variations in mass density, pressure, 
and stellar radius, becomes localized inside the star. In the subsequent discussion, we will show that shifts in QNM 
frequencies can be expressed in terms of integrals with finite domains of integration. 



B. First order perturbation theory 



Consider an unperturbed neutron star which has a scaled circumferential (tortoise) radius Rq (R*q), a scaled 
potential Vo(r»), and a scaled known mass distribution function mo(r). Note that all physical quantities are scaled 
with its own mass Mq. The unperturbed star (e.g. TVIIM) will be used to provide the zeroth-order approximation 
of a realistic star whose physical parameters are not exactly known and to be determined from the inversion scheme. 

Analogous to other conventional perturbation theories, we will assume that the QNM wave function of the unper- 
turbed star is known. If now the star is perturbed to a new configuration which is identical to that of the realistic 
star in consideration, the new potential function is given by: 



(5.14) 



As a consequence, the circumferential (tortoise) radius of the perturbed star is also modified such that, R = Rq + R\ 
(i?* = i?*o + i?*i). However, outside the star, V, as a function of f, remains unchanged. 

Consider the q-th QNM of the unperturbed star, which has a scaled frequency uj q °\ From SCLPT we obtain the 
first-order shift in the scaled frequency ui = Mm, which reads |35j: 



(5.15) 



with 



R*i 



d/, 



(0) 



df* 



Ri 



df 



(5.16) 



and i/> g (r*) being the unperturbed wave function of the q-th QNM mode. Besides, fq°\r*) is the logarithmic derivative 
of the unperturbed wave function ipqi?*) that satisfies both the regular boundary condition at f* = and the outgoing 

wave boundary condition at spatial infinity, while fq + \r*, u>) is the logarithmic derivative of the outgoing wave solution 
of frequency ui to the scaled NSRW equation 33]. From Q5.15[) . (|5.1t)|) and (|5.17|) it is then obvious that one has to 
evaluate the changes in the scaled potential, the scaled circumferential and tortoise radii of the star in order to find 

Uq . 
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C. The perturbed star 

The perturbed star is in fact the realistic star considered in the inversion scheme, whose scaled mass distribution 
m(f) and density distribution /5(f) are expressed in terms of f and its own exact compactness C = M/R as follows: 



m(f) 



m(r] 
M 



p-2 

m (f,C) + ^2fi j m J (r,C), 

3 = 1 

p-2 

p(r) = M 2 p(r) =p (f,C) + J2fi j p j (r,C). 



(5.18) 
(5.19) 



Here 



m Q (f,C) 



mo(r,M,R) 1 , „, ,, „r r, 
M = f - 3C ? )> 

15 



p (r,C) = M 2 p (r,M,R) = -fC 3 (l-C 2 ? 2 ), 

&(f,C) = M 2 p 3 (r, M, R) = ^ +1 ^- 2 (C 2 r 2 - 1) [(5 + 2j)C 2 f 2 - (1 + 2j)] . 

In order to compare the perturbed and unperturbed stars, we rewrite fh{f) and p{f) as: 

m(f) - m (0) (f)+m (1) (f), 
/5(r) - /5(°)(f)+p( 1 )(f), 



with 



m (0) (f) - mo(f,Co), 

p<°>(f) = p (r,Co), 
p-i 

m (1) (f) = ^/ij-mj^Co), 
p-i 

p«(f) = ^^(f.Co), 
j'=i 



(5.20) 
(5.21) 
(5.22) 
(5.23) 



(5.24) 
(5.25) 



(5.26) 
(5.27) 

(5.28) 
(5.29) 



= y (C 2 r 3 -C 4 f 5 ), 



p P -i(0 



C=Cn 
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8^ 



2^2 \ 



Co(3-5^r~ 



(5.30) 
(5.31) 



and /Up-i = C — Co = C^ 1 ' as mentioned previously. m^(f) and p^^f) are considered as the first order changes in 
the mass and density distributions. Besides, it is straightforward to show that 



R Rq 



■ C (o)> 



(5.32) 



R^=R^exp 



i>W(RW)-\W(RW) 



+ 



Jo 



exp 



A(°)(f)-^°)(f) 



dr. 



(5.33) 
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These first-order quantities, including to' 1 ), R* , are the inputs to SCLPT. Once ujq 1 ^ is found from the 

first-order SCLPT, the first order change in QNM frequency, a/ 1 ), can be obtained from the following relation: 



,~.(0) , -(1) - (0) 

(1) uJq T LOq UJq 



M(°) + M« M(°) 

d 1} mVq^ 



Af(o) [M(°)f 



(5.34) 



D. Perturbed TOV equation 

The perturbation in m^(f) (or p^(r)) of course leads to variations in pressure and metric coefficients, and in turn 
affects the potential V. So we need to consider the scaled TOV equations [s^lssT], 



dm 
df 



47rf 2 p(f), 



dv _ 2m(f) + 8irr 3 P(r) 
df f(f — 2m(f)) ' 

dP 1 . . , ~,ssdv 

with v = v . Denoting the first order changes in P and v by P^ and we can show from the scaled TOV equations 
that: 



df 
df 



1 



(f- 2m(°)) 2 
1 



(2 + 167rf 2 p(°))m (1) + 8?rf 2 (f - 2#')pW 
dvW , n , -, nV ^W 



(p«+pw)^ r + (p(°) + p( ))^ 



(5.35) 
(5.36) 



Therefore, both P^ and can be evaluated by numerically solving this pair of differential equations. Besides, 
from 



-A(f) 



1 



2m(f) 



(5.37) 



we have 



_ 2m( 1 )exp(AW) 



(5.38) 



E. Perturbed potential 

Upon the changes in mass and density functions, the scaled V rw varies in accordance with (I5.14|) . However, V is 
an explicit function of f, instead of f*. In order to express V as a function of f* such that our perturbation scheme is 
applicable, we have to regard f as a function of f* and in turn obtain the implicit function V(f(f*)). In the absence 
of perturbation, f = f(°)(f„) is determined from the defining equation of the tortoise coordinate: 

exp{[-^(f) + \^(f)]/2}df. (5.39) 

With the introduction of fhS x > and the scaled circumferential radius is modified such that f = f^'if*) + f( x \f*), 
where f^(f^) is defined via the following equation: 

r( '(f»)+f (1) (P») 

exp{[A (0) (f) + A (1) (f) - £ (0) (f) - u {1 \f)]/2}df. (5.40) 
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Neglecting second order terms, we can find ' in terms of f*: 



r^ 1 ' (f*) = exp 



By virtue of (|5.14|1 . we have 



>(?») 



exp 



dr. 



(5.41) 



7(f(°))-y(°)(f( )) 



dm ' 



(5.42) 



and from (|5.12|) explicit expressions for the terms in the right hand side of (|5.42|1 can be found: 



dVW f r(0 )di>W p m f_ 21(1 + 1) / rfp (0) dP (0) \ 6[47r(f W ) 3 /5(°) - 3m' )] | 

<ff(°) ~ 9f(°) +C | (f(°)) 3 + ^dfW df(°) ) (r(°)) 4 J' 



(5.43) 



V"(f(o))_v-(o)(f(o)) = jj(i)7(o)(f(o)) + e *<' 



47T 



(f(0))3 



(5.44) 



The results obtained above together with the first order SCLPT lead explicit expressions for c q j in (|4.1|) 



VI. NUMERICAL RESULTS 



In Figs. 151- 1101 we show numerical results obtained from the inversion scheme developed in the present paper for six 
different realistic neutron stars with a common compactness of 0.28. To gauge the accuracy of the scheme, we have 
used four different combinations of QNM frequencies: (i) w\\ (ii) u>\ and ui2\ (hi) oj\, 102 and W3; and (iv) u)\ and uj\\. 
Here 0^1,^2,^3, ■ ■ ■ are ordered in increasing frequencies, and uiu is the frequency of a wn mode |40l |. In all cases the 
number of parameters p was chosen to be twice the number of QNMs (i.e. In). It is clearly seen from these figures 
that (i) the results obtained from using only the frequency of the least damped mode, u>\ (empty symbols), can nicely 
reproduce the EOS except at the high density regime; (ii) the results obtained from a combination of u>\ and L02 (grey 
symbols) are readily improved; (iii) using all three frequencies w\, LO2 and W3 (dark symbols) indeed yields a perfect 
match with the exact values; and (iv) the results obtained from the combination of u>i and u>u (crosses) are slightly 
less accurate than those from the combination of lo\ and ll>2. 

In Fig. 13 we compare these six EOSs and demonstrate that it suffices to use only two frequencies {ui\ and LO2) 
to accurately reproduce and distinguish these EOSs. In Table I we tabulate the exact values of M and R and 
those obtained from our scheme using one or two QNMs for stars constructed with these EOSs. The result is truly 
encouraging, especially for the mass, which is very close to the exact value in all cases. Despite that the computed 
radius might deviate from the exact one by a few percent in some cases, the result is in fact very accurate taking into 
account of the fact that the mass density near the surface of realistic neutron stars is usually very low. Hence, we 
conclude that the frequencies of two leading w-mode oscillations can readily lead to accurate determination of M, R 
and EOS. 

In each inversion process the number of parameters p is less than or equal to twice the number of QNMs used (i.e. 
2ra). This can be easily understood as each QNM frequency carries two data, namely the pulsation and decay rates. 
As a result, at most 2n parameters can be determined from the inversion process. On the other hand, as shown in 
Fig. 1111 where the inversion scheme is applied to an APR1 star with three known QNMs, the case p = 6 obviously 
yields the best result. Therefore, for a fixed number of QNMs, the choice p = 2n seems to be the optimal one and 
this remark is supported by our numerical results. 

The stability of the inversion process has been investigated by introducing artificial noises into the real and imaginary 
parts of QNM frequencies. In each inversion process the radius of the star in consideration, which was also an unknown 
to be determined from the process, was divided into a fixed number of equal partitions. The values of the radius, 
density, and pressure at each of these small partitions were in turn obtained from the inversion scheme. After gathering 
the data, we averaged them over a number of simulations in which random noises were superimposed on the values 
of QNM frequencies, and evaluated the respective error bars. Fig. ^| shows the result for an APR1 star with two 
approximately known QNM frequencies, lo\ and 0J2, and four parameters were used in the inversion. The random 
noise introduced was within ±5% of the exact value. The grey circles (with error bars) denote the average values 
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of 12 simulations, which are fairly accurate. This strongly suggests that the inversion scheme proposed here is still 
feasible even if gravitational signals are blurred by experimental noises. However, there was a need for caution to be 
exercised in inverting noise-polluted QNM frequencies. In some cases the result obtained might explicitly violate the 
thermodynamic stability condition dp /dp. We have in fact discarded such ill-defined data in the inversion scheme 
mentioned above. 

Lastly, to further demonstrate the robustness of the inversion scheme, we have also applied our scheme to infer the 
EOS of a SQS @. The EOS of quark matter is given by the MIT bag model \ri& : 

P=(p-4B)/3, (6.1) 

where the value of the bag constant B is 57MeVfm -3 . The crust of the SQS star, a thin layer of normal matter 
obeying the Baym-Pethick-Sutherland EOS is supported by Coulomb forces above the quark matter core [l3j . 
The EOS for quark matter is rather exotic and the density profile of such a star is even discontinuous across the inner 
boundary of the crust. However, as shown in Fig. 1131 the inversion scheme still works nicely as long as four QNMs 
are used. It is then obvious that the current inversion scheme outperform the one proposed by Lindblom, which is 
limited to barotropic type EOS [34| . 



VII. DETECTION OF QNMS 



In the above discussion we show how the internal structure of a neutron star can be systematically determined 
from several of its QNM frequencies. Like all other studies in the area of gravitational-wave asteroseismology (see 
e -g- [2ll I22I I231 I24I EM l26l I27I l28l |29|). our method depends crucially on reliable data of gravitational wave signals, 
which have not yet been available. Therefore, the realization of the present scheme hinges on the successful operation 
of various existing gravitational-wave detectors (such as LIGO, VIRGO, GEO600), and also the development in the 
design and construction of new gravitational-wave detectors in the future. In the following discussion, we consider 
the possibility of detecting QNMs of neutron stars based on analysis outlined in Ref. [24| • 

The detection of a QNM signal, an exponentially decaying sinusoidal wave with frequency / = oj r /(27r) and decay 
time r = 1/oJi, depends on the distance r between the emitter (the neutron star) and the observer (the earth), the 
energy E gw channelled in such a QNM, and of course the sensitivity of the gravitational-wave detector, which is 
measured in terms of its noise power spectral density S n (/) (see Fig. I14|l . A strong signal is associated with a large 
signal-to-noise ratio given by [24( : 



J l + 4Q 2 2S n ' 

where the quality factor Q = 7t/t, and the amplitude of such a gravitational wave signal |24| : 

£ gw \ 1/2 ^10kpc\ ^lkHz^ ^lms x 1/2 



^"-°-"(twJ {-^){— ){-) ■ <"> 

The fractional statistical errors in the determination of the frequency and the decay time, respectively denoted by 
<7f/f and oy/t, are estimated to be [24j : 



07 , /l-2Q 2 + 8Q 4 

' ~ 0.0042P- 1 ' ^ 



/ V 4 Q 4 



(^-r\ (7.3) 

V 1 ms / 



--0.013P-S (7.4) 



with 



t \ Q 2 VlkHz 



E \~ 1/2 / r \ ( <? 1/2 



P Uo- 6 M Q c 2 J UokpcJ (^0-23 Hz" 1 / 2 J ' {7 ' 5) 

Both LIGO and VIRGO have a broad-band sensitivity in the hundred-Hertz range where Sl /2 - 1(T 23 Hz~ 1/2 
(see Fig. 1140 . It is then obvious that these first generation detectors could not accurately locate ui-mode QNMs 



12 



of neutron stars in our galaxy, which are characterized by the following typical values: / ~ 10 kHz, r ~ 0.02 ms, 
Egw ~ 10 _6 M Q c 2 and r ~ lOkpc (24|. In fact, it is generally believed that the sensitivities of the first generation 
detectors are not sufficient to carry out quantitative measurements. 

LIGO-II, a second generation gravitational-wave detector, is aimed at operation in the hundred-Hertz range near 
the standard quantum limit (SQL) and has a lower noise level given approximately by |4l| : 



Sn(f) 

Hz- 1 



= 10 



-49 



x- 4 - 14 - 5x~ 2 • 



(111)(2- 2x 2 + x 4 ) 
Y+x~ 2 



(7.6) 



with x = //(215 Hz). In spite of the improved sensitivity, as shown in Table ITU where E gw = 1O _6 M0C 2 and 
r = 10 kpc are assumed, LIGO-II is not yet sensitive enough to locate w-mode QNMs of neutron stars in our own 
galaxy. On the other hand, it is worthwhile to note that SQL can in fact be overcome by utilizing quantum non- 
demolition (QND) devices [42]. The design of LIGO-III, a third generation gravitational-wave detector, is currently 
under intense discussion. Various strategies that can beat the gravitational-wave SQL by an arbitrarily large amount, 
over an arbitrarily wide range of frequencies have been proposed and analyzed |43l |44| . 

In the present paper we consider specifically another third generation gravitational-wave detector, EURO (European 
Gravitational Wave Observatory), proposed by a group of European scientists in 2000, aims to make quantitative 
surveys of gravitational- wave signals over 4 decades of frequency, namely from 1 Hz to tens of kHz [29t Elj . As 
stated in the proposal |4l| |. "EURO can do asteroseismology on neutron stars by resolving the frequencies of the 
dominant normal modes of vibration". In the proposal two possible configurations of EURO were discussed. In the 
first configuration, EURO, the sensitivity in the high frequency range is still limited by the photon shot noises, with 
a noise level given by [4l| : 



Sn(f) 

Hz- 1 



10 



-50 



3.6 x 10 9 1.3 x 10 5 



f 4 



r- 



!■•'! • ID" 7/. (1 + 4 

Jk 



(7.7) 



where the knee-frequency fj~ = 1000 Hz. As shown in Fig. O the noise level of this configuration in the kHz range 
is two orders less than that of LIGO-II. With this detector w-mode gravitational wave signals emitted from neutron 
stars within our galaxy can be detected with error bound of a few percent. Therefore, the inversion scheme outlined 
above will become viable once EURO operates as proposed. 

In the second configuration, referred to as EURO (xylophone) in the present paper, the shot noise limit in the 
high-frequency regime is beaten by parallel operation of several narrow-banded cryogenic interferometers comprising 
a "gravitational- wave xylophone" [23, El • The noise level of such a device is approximately given by ^3] : 



Sn(f) 

Hz- 1 



10 



-50 



3.6 x 10 9 1.3 x 10 5 



f 4 



P 



(7.8) 



which can be obtained by removing the last term (i.e. the shot-noise term) in (17. 7|) . As can be inferred from Fig. 1141 
Table [HJ and also the following scaling relations |2j 




Ea 



(7.9) 
(7.10) 
(7.11) 



the advent of the detector EURO (xylophone) will enable observation of gravitational waves from extra-galactic 
neutron stars and black holes. We therefore expect that the inversion scheme developed in the present paper can reveal 
fruitful information about the internal structure of neutron stars from the data gathered by these third generation 
gravitational-wave detectors. 



VIII. CONCLUSION 



Despite the mass of some neutron star binaries can be inferred from their orbital periods, as yet there is no generic 
method to determine their radii except for a few cases where data of gravitational red-shift are available 0, ■ In the 
present paper we have proposed a robust inversion scheme to determine the mass, the radius, the mass distribution 
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and the EOS of a neutron star from its gravitational wave spectra, which could be obtained by the matched filtering 
method (see, e.g. HU). We expect that our scheme could operate in conjunction with gravitational- wave detectors 
available in one or two decades, e.g. LIGO-III, EURO and EURO (xylophone), to probe the interior of neutron stars. 
As clearly shown in Fig. [21 an inversion scheme using two leading QNM frequencies can readily discriminate between 
any two smooth EOSs considered there. The method developed in the present paper in fact makes the construction 
of such gravitational-wave detectors more rewarding. 

In reality, gravitational signals are inevitably contaminated by noises of various kinds, including thermal, quantum 
and gravity-gradient noises (see, e.g. 0> HE HE EH HE E3 an( A references therein). Analogous to other fields of 
astronomy, observational errors introduced by noises are deemed the major hurdle to the implementation of the inver- 
sion scheme proposed here. In order to infer the fine internal structure of neutron stars, theoretical and experimental 
efforts must be done to minimize the noise level, and sophisticated data acquisition skills have to be developed to 
distinguish signals from noises [ill \I§L |44|. On the other hand, as demonstrated in Fig. ^1 our method remains 
accurate and stable against contamination of data created by noises. The construction and operation of advanced 
gravitational-wave detectors that can provide higher detection sensitivities and resolution power will strongly boost 
the progress in gravitational-wave asteroseismology [SEE El The proposal in the current paper will then become 
a powerful tool. As predicted in Refs. [23,|4l|, the prospect of gravitational-wave asteroseismology is bright. 
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Ivl 


IX I Kill J 


APRl 


4.460/4.475/4.460 


11.83/12.12/11.70 


APR2 


4.173/4.184/4.173 


11.07/11.26/10.91 


AU 


3.843/3.862/3.843 


10.19/10.61/10.08 


A 


3.289/3.282/3.289 


8.723/8.529/8.564 


C 


2.791/2.788/2.791 


7.404/7.185/7.216 


UT 


3.654/3.651/3.655 


9.691/9.745/9.664 



TABLE I: The mass M and the radius R of neutron stars constructed with different EOSs are compared with the corresponding 
values obtained from inversion using one or two QNMs. In each entry the values are obtained from TOV equation/inversion 
using uji /inversion using w\ and u)i respectively. 



Detector 


/ (kHz) 


r (ms) 


[S„(/)] 1/2 (Hz- 1/2 ) 


(syAO 2 


°f/f 


<J t /t 


LIGO II 


7.320 
12.76 


0.065 
0.037 


1.13 x 10"^ 
1.98 x 10" 22 


3.76 x 10" 1 
4.06 x 10" 2 


0.99 
3.07 


3.8 
11.6 


EURO 


7.320 
12.76 


0.065 
0.037 


8.42 x 10" 2b 
1.46 x 10" 24 


6.80 x 10 d 
7.45 x 10 2 


7.4 x 10" J 
2.3 x 10" 2 


2.8 x 10'^ 
8.5 x 10" 2 


EURO 
(xylophone) 


7.320 
12.76 


0.065 
0.037 


4.93 x lO" 2 ' 
2.83 x 10" 27 


1.99 x 10 K 
1.99 x 10 8 


4.3 x 10"" 

4.4 x 10" 5 


1.7 x 10"* 
1.7 x 10" 4 



TABLE II: The noise amplitude spectral density [Sn(f)] 1 ^ 2 (Hz~ 1//2 ), the signal to noise ratio (S/N) 2 , o//f and ov/r are 
tabulated for LIGO II, EURO and EURO (xylophone) gravitational-wave detectors. Here the frequency and damping time of 
the two leading axial w-modes of an APRl star with C = 0.28 are considered, which are also used as the input in Fig. 1121 and 
£ gw = 10" 6 M Q c 2 and r = lOkpc are assumed. 



16 



0.50 



CD 
DC 




0.20 

M/R 



0.30 



FIG. 1: The real and imaginary parts of Mu> for the least-damped axial w-mode of six realistic stars (unfilled symbols) and 
TVIIM (stars) are shown as a function of M/R in panels (a) and (b) respectively. The dotted line represents the best quadratic 
fit to those of TVIIM. 




FIG. 2: Six different EOSs are shown by lines of different formats. Grey symbols are the corresponding values obtained from 
inversion scheme using uii and ui2 for stars with C = 0.28. Here P and p are in cgs units. 
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FIG. 3: The scaled mass distributions m(r)/M of various realistic neutron stars and TVIIM with a common compactness 
C = 0.28 are plotted against r/R. 
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FIG. 4: Re(w (c) )/Im(w (c) ) for the three leading modes of TVIIM star (i) u[ c > (solid-line); (ii) lo { 2 c) (dashed-line) and (iii) u 
(dot-dashed-line) are plotted against the compactness M/R. 
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FIG. 5: Panels (a) and (b) depict p(r) and P(p) respectively for an APR1 star with C = 0.28. The solid line is the theoretical 
value, while the unfilled/grey /dark circles represent the results obtained from inversion scheme using one/two/three leading 
axial w-modes. The result obtained from inversion scheme using u)\ and the frequency of a lun-mode is shown by the crosses. 
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FIG. 6: Same as Fig. EI except that an APR2 star is considered instead. 
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FIG. 7: Same as Fig. Q3 except that an AU star is considered instead. 




FIG. 8: Same as Fig. except that a model A star is considered instead. 
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FIG. 9: Same as Fig.0 except that a model C star is considered instead. 
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FIG. 10: Same as Fig. except that a UT star is considered instead. 
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FIG. 11: Panels (a) and (b) depict p(r) and P(p) respectively for an APR1 star with C = 0.28. The solid line is the theoretical 
value, while the left triangles/squares/dark circles represent the results obtained from inversion scheme using four/five/six 
parameters and three leading axial w-modes. 
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FIG. 12: Panels (a) and (b) depict p(r) and P(p) respectively for an APR1 star with C = 0.28. The solid line is the 
theoretical value, while the grey circles connected by a grey line represent the results obtained from inversion scheme using two 
approximately known leading axial ui-modes and four parameters. Twelve sets of data were used in the simulation. The errors 
in the real and imaginary parts of the QNM frequency were uniformly distributed in a ±5% interval. 
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FIG. 13: Panels (a) and (b) depict p(r) and P(p) respectively for a SQS star with C = 0.2. The solid line is the theoretical 
value, while the grey triangles/dark triangles/crosses respectively represent the results obtained from inversion scheme using 
two, three and four leading axial w-modes. 
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FIG. 14: The noise amplitude Sn (/) (Hz -1 ^ 2 ) is plotted against / (Hz) for gravitational- wave detectors LIGO-I, LIGO- 
II, EURO and EURO (xylophone). The curves are based on the interpolation formulas available at the EURO homepage 
http://www.astro.cf.ac.uk/geo/euro 



